Chapter 5 HMSC
5.1 Setup
# Random effects data (study design)
StudyDesign <- sample_metadata %>%
select(sample,Location) %>%
rename(location=Location) %>%
column_to_rownames("sample") %>%
mutate(location = factor(location))
# Genome count table (quantitative community data)
YData <- read_counts %>%
mutate(across(where(is.numeric), ~ . +1 )) %>% #add +1 pseudocount to remove zeros
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / 150) )) %>% #transform to genome counts
mutate(across(where(is.numeric), ~ log(.) )) %>% #log-transform
arrange(genome) %>%
column_to_rownames("genome") %>%
select(all_of(row.names(StudyDesign))) %>% #filter only faecal samples
as.data.frame() %>%
t() # transpose
# Fixed effects data (explanatory variables)
XData <- sample_metadata %>%
column_to_rownames("sample") %>%
mutate(logseqdepth=read_counts %>% #total log-sequencing depth
select(all_of(row.names(StudyDesign))) %>%
colSums() %>%
log()
) %>%
rename(origin=Origin, sex=F.M) %>%
mutate(origin = factor(origin)) %>%
mutate(sex = factor(sex)) %>%
select(origin, sex, logseqdepth)
# Genome phylogeny
PData <- genome_tree5.3 Define and Hmsc models
#Define models
model1 = Hmsc(Y=YData,
XData = XData,
XFormula = XFormula,
studyDesign = StudyDesign,
phyloTree = PData,
ranLevels=list("location"=rL.location),
distr = "normal",
YScale = TRUE)
#Save list of models as an R object.
model_list = list(model1=model1)
if (!dir.exists("hmsc")){dir.create("hmsc")}
save(model_list, file = "hmsc/hmsc.Rdata")Upload hmsc/hmsc.Rdata to the HPC respecting the directory structure.
5.5 Generate Hmsc executables
The next chunk generates shell files for every combination of model, MCMC samples and MCMM thinning, ready to be launched as SLURM jobs.
modelchains <- expand.grid(model = names(model_list), sample = MCMC_samples_list, thin = MCMC_thin_list)
if (!dir.exists("hmsc")){dir.create("hmsc")}
for(i in c(1:nrow(modelchains))){
modelname=as.character(modelchains[i,1])
sample=modelchains[i,2]
thin=modelchains[i,3]
executablename <- paste0("hmsc/exe_",modelname,"_",sample,"_",thin,".sh")
fitname <- paste0("fit_",modelname,"_",sample,"_",thin,".Rdata")
convname <- paste0("conv_",modelname,"_",sample,"_",thin,".Rdata")
model <- paste0('model_list$',modelname)
psrf.beta.name <- paste0("psrf.beta.",modelname,"_",sample,"_",thin)
psrf.gamma.name <- paste0("psrf.gamma.",modelname,"_",sample,"_",thin)
psrf.rho.name <- paste0("psrf.rho.",modelname,"_",sample,"_",thin)
jobname <- paste0("hmsc_",modelname,"_",sample,"_",thin)
minutes <- 1000
code <- sprintf("#!/bin/bash
#SBATCH --job-name=%s # Job name
#SBATCH --nodes=1
#SBATCH --ntasks=4 # Run on 4 CPUs
#SBATCH --mail-user=antton.alberdi@sund.ku.dk
#SBATCH --mem=96gb # Job memory request
#SBATCH --time=%d # In minutes
# Activate conda environment
module load mamba/1.3.1
source activate /maps/projects/mjolnir1/people/jpl786/AMAC001_fibre_trial/hmsc/hmsc_env
# Run R script
Rscript -e '
library(tidyverse)
library(Hmsc)
# Load formulas and data
load(\"hmsc.Rdata\")
# Declare placeholders
modelname = \"%s\"
model = %s
fitname = \"%s\"
convname = \"%s\"
sample = %d
thin = %d
nchains = %d
# Run model fitting
m = sampleMcmc(hM = model,
samples = sample,
thin = thin,
adaptNf=rep(ceiling(0.4*sample*thin),model$nr),
transient = ceiling(0.5*sample*thin),
nChains = nchains,
nParallel = nchains)
# Assess chain convergence
mpost = convertToCodaObject(m,
spNamesNumbers = c(T,F),
covNamesNumbers = c(T,F),
Beta = TRUE,
Gamma = TRUE,
V = FALSE,
Sigma = FALSE,
Rho = TRUE,
Eta = FALSE,
Lambda = FALSE,
Alpha = FALSE,
Omega = FALSE,
Psi = FALSE,
Delta = FALSE) # Convert to CODA object
# Fixed effects
assign(paste0(\"psrf.beta.\", modelname,\"_\",sample,\"_\",thin), gelman.diag(mpost$Beta,multivariate=FALSE)$psrf)
# Traits
assign(paste0(\"psrf.gamma.\", modelname,\"_\",sample,\"_\",thin), gelman.diag(mpost$Gamma,multivariate=FALSE)$psrf)
# Phylogeny
assign(paste0(\"psrf.rho.\", modelname,\"_\",sample,\"_\",thin), gelman.diag(mpost$Rho,multivariate=FALSE)$psrf)
# Write convergence data
save(%s, %s, %s, file=convname)
# Save model fit object
save(m, file=fitname)
'
", jobname, minutes, modelname, model, fitname, convname, sample, thin, nChains, psrf.beta.name, psrf.gamma.name, psrf.rho.name)
writeLines(code, executablename)
}Upload the produced hmsc/exe_XXXXX.sh files to the HPC respecting the directory structure.
5.8 Variance partitioning
# Compute variance partitioning
varpart=computeVariancePartitioning(m)
varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
group_by(variable) %>%
summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
tt()| variable | mean | sd |
|---|---|---|
| Random: location | 37.900015 | 25.317903 |
| logseqdepth | 56.110626 | 25.796874 |
| sex | 4.937460 | 5.612719 |
| origin | 1.051899 | 1.282563 |
# Basal tree
varpart_tree <- genome_tree
#Varpart table
varpart_table <- varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))
#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% varpart_tree$tip.label) %>%
arrange(match(genome, varpart_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% varpart_tree$tip.label) %>%
arrange(match(genome, varpart_tree$tip.label)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
select(colors) %>%
pull()
# Basal ggtree
varpart_tree <- varpart_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()
# Add variance stacked barplot
vertical_tree <- varpart_tree +
scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
geom_fruit(
data=varpart_table,
geom=geom_bar,
mapping = aes(x=value, y=genome, fill=variable, group=variable),
pwidth = 2,
offset = 0.05,
width= 1,
orientation="y",
stat="identity")+
labs(fill="Variable")
vertical_tree5.9 Posterior estimates
# Select desired support threshold
support=0.9
negsupport=1-support
# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
mutate(value = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
pivot_wider(names_from = variable, values_from = value) %>%
#select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
column_to_rownames(var="genome")
#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
select(colors) %>%
pull()
# Basal ggtree
postestimates_tree <- postestimates_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()
# Add posterior significant heatmap
postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
labs(fill="Trend")
postestimates_tree +
vexpand(.25, 1) # expand top
## Correlations
#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)
# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
keep.tip(., tip=m$spNames)
#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
+ (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean
matrix <- toPlot %>%
as.data.frame() %>%
rownames_to_column(var="genome1") %>%
pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
ggplot(aes(x = genome1, y = genome2, fill = cor)) +
geom_tile() +
scale_fill_gradient2(low = "#be3e2b",
mid = "#f4f4f4",
high = "#b2b530")+
theme_void()
htree <- genome_tree_subset %>%
force.ultrametric(.,method="extend") %>%
ggtree(.)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#create composite figure
grid.arrange(grobs = list(matrix,vtree),
layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1)))5.10 Predict responses
# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")
gradient = c("domestic","feral")
gradientlength = length(gradient)
#Treatment-specific gradient predictions
pred <- constructGradient(m,
focalVariable = "origin",
non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
predict(m, Gradient = ., expected = TRUE) %>%
do.call(rbind,.) %>%
as.data.frame() %>%
mutate(origin=rep(gradient,1000)) %>%
pivot_longer(!origin,names_to = "genome", values_to = "value")# weights: 9 (4 variable)
initial value 101.072331
final value 91.392443
converged
5.10.0.1 Function level
functions_table <- elements_table %>%
to.functions(., GIFT_db) %>%
as.data.frame()
community_functions <- pred %>%
group_by(origin, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
select(-row_id) %>%
column_to_rownames(var = "origin") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(functions_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="origin")
})#max-min option
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
function_predictions <- map_dfc(community_functions, function(mat) {
mat %>%
column_to_rownames(var = "origin") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_functions[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)
# Positively associated
function_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D02 | 8.355220e-03 | -0.0037113857 | 0.0205920876 | 0.818 | 0.182 |
| B08 | 7.253915e-03 | -0.0038460307 | 0.0187179563 | 0.785 | 0.215 |
| B01 | 1.006611e-03 | -0.0062300881 | 0.0080852143 | 0.603 | 0.397 |
| S01 | 9.339371e-04 | -0.0117178266 | 0.0135109195 | 0.591 | 0.409 |
| B09 | 1.304867e-05 | -0.0005590059 | 0.0004391979 | 0.356 | 0.644 |
# Negatively associated
function_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D08 | -0.0011620023 | -0.002196111 | -0.0002312023 | 0.041 | 0.959 |
| B03 | -0.0108172543 | -0.018252135 | -0.0035018570 | 0.046 | 0.954 |
| D06 | -0.0032232182 | -0.006932327 | 0.0001148282 | 0.104 | 0.896 |
| B04 | -0.0083179762 | -0.018814251 | 0.0009380017 | 0.131 | 0.869 |
| D07 | -0.0117417980 | -0.029291626 | 0.0042684819 | 0.184 | 0.816 |
| B06 | -0.0064687006 | -0.016619760 | 0.0033329743 | 0.196 | 0.804 |
| D05 | -0.0019072483 | -0.007708018 | 0.0032904680 | 0.218 | 0.782 |
| D03 | -0.0041396136 | -0.013006620 | 0.0033421571 | 0.230 | 0.770 |
| S03 | -0.0097522350 | -0.033201541 | 0.0148288014 | 0.248 | 0.752 |
| B02 | -0.0037333239 | -0.012806201 | 0.0054099794 | 0.283 | 0.717 |
| D09 | -0.0020780433 | -0.008498140 | 0.0048655860 | 0.312 | 0.688 |
| S02 | -0.0043222701 | -0.013090770 | 0.0034865384 | 0.323 | 0.677 |
| B07 | -0.0039506685 | -0.015755996 | 0.0084548853 | 0.342 | 0.658 |
| D01 | -0.0003552259 | -0.005071442 | 0.0042893404 | 0.411 | 0.589 |
| B10 | -0.0000145377 | -0.000309013 | 0.0002678120 | 0.477 | 0.523 |
positive <- function_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
select(-negative_support) %>%
rename(support=positive_support)
negative <- function_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.02,0.02)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")5.10.0.2 Element level
elements_table <- genome_gifts %>%
to.elements(., GIFT_db) %>%
as.data.frame()
community_elements <- pred %>%
group_by(origin, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
select(-row_id) %>%
column_to_rownames(var = "origin") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(elements_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="origin")
})
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
element_predictions <- map_dfc(community_elements, function(mat) {
mat %>%
column_to_rownames(var = "origin") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_elements[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)# Positively associated
element_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
filter(positive_support>=0.9) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D0205 | 0.012685988 | 0.0026218528 | 0.022934401 | 0.956 | 0.044 |
| D0208 | 0.010171843 | 0.0016293321 | 0.018293536 | 0.941 | 0.059 |
| D0906 | 0.003664187 | 0.0001651717 | 0.007903952 | 0.929 | 0.071 |
| B0103 | 0.008392278 | 0.0009118417 | 0.016138328 | 0.924 | 0.076 |
| D0507 | 0.003661915 | 0.0001780608 | 0.007286307 | 0.910 | 0.090 |
| D0504 | 0.004349264 | 0.0000896653 | 0.009239292 | 0.908 | 0.092 |
element_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
filter(negative_support>=0.9) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D0801 | -0.001648434 | -0.002310322 | -9.380558e-05 | 0.013 | 0.987 |
| D0802 | -0.001648434 | -0.002310322 | -9.380558e-05 | 0.013 | 0.987 |
| B0709 | -0.002156801 | -0.003595694 | -6.174761e-04 | 0.028 | 0.972 |
| D0517 | -0.004675577 | -0.007725113 | -1.322827e-03 | 0.028 | 0.972 |
| B0302 | -0.005012428 | -0.010123378 | -5.918431e-04 | 0.029 | 0.971 |
| D0603 | -0.001993397 | -0.003715990 | -4.415768e-04 | 0.032 | 0.968 |
| D0601 | -0.009415187 | -0.016646322 | -2.368741e-03 | 0.033 | 0.967 |
| B0310 | -0.013036483 | -0.024177638 | -3.420713e-03 | 0.038 | 0.962 |
| D0611 | -0.004154585 | -0.009328677 | -2.659067e-04 | 0.039 | 0.961 |
| D0903 | -0.004154585 | -0.009328677 | -2.659067e-04 | 0.039 | 0.961 |
| B0219 | -0.004197892 | -0.009466696 | -2.638587e-04 | 0.040 | 0.960 |
| D0610 | -0.003089315 | -0.004942868 | -8.590712e-04 | 0.043 | 0.957 |
| D0817 | -0.005054593 | -0.010663544 | -5.732221e-04 | 0.052 | 0.948 |
| D0606 | -0.005946253 | -0.010886608 | -1.062497e-03 | 0.053 | 0.947 |
| D0807 | -0.004229573 | -0.008646531 | -7.602551e-04 | 0.055 | 0.945 |
| B0214 | -0.021448686 | -0.040019414 | -4.451973e-03 | 0.060 | 0.940 |
| D0908 | -0.015981409 | -0.027608998 | -3.342110e-03 | 0.062 | 0.938 |
| B0303 | -0.011717078 | -0.021165665 | -2.389863e-03 | 0.063 | 0.937 |
| B0804 | -0.016699845 | -0.031709184 | -2.626835e-03 | 0.066 | 0.934 |
| B0601 | -0.009487153 | -0.018010058 | -1.308593e-03 | 0.071 | 0.929 |
| D0612 | -0.001696748 | -0.003105924 | -2.115429e-04 | 0.072 | 0.928 |
| B0603 | -0.016221363 | -0.030525848 | -1.797226e-03 | 0.073 | 0.927 |
| B0401 | -0.011912922 | -0.024386151 | -7.694553e-04 | 0.076 | 0.924 |
| D0508 | -0.003348671 | -0.007220707 | -1.785349e-04 | 0.082 | 0.918 |
| D0816 | -0.005917943 | -0.012001822 | -3.994812e-04 | 0.084 | 0.916 |
| D0512 | -0.018031053 | -0.037115248 | -1.288607e-03 | 0.087 | 0.913 |
| B0309 | -0.007714059 | -0.015654662 | -3.596295e-04 | 0.094 | 0.906 |
| B0708 | -0.024244082 | -0.047246442 | -7.506292e-06 | 0.100 | 0.900 |
positive <- element_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
select(-negative_support) %>%
rename(support=positive_support)
negative <- element_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.04,0.04)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")